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Introduction 

Considerable progress over the past thirty years has 
been made in the development, of large-scale computa- 
tional fluid dynamic s (CFD) solvers for the Euler and 
Navier- Stokes equations. Computations are used rou- 
tinely to design the cruise shapes of transport aircraft 
through complex-geometry simulations involving the 
solution of 25-100 million equations; in this arena, the 
number of wind-tunnel tests for a new design has been 
substantially reduced. However, simulations of the en- 
tire flight envelope of the vehicle 1 , including maximum 
lift, buffet onset, flutter, and control effectiveness, have 
not been as successful in eliminating the reliance on 
wind-tunnel testing. These simulations involve un- 
steady flows with more separation and stronger shock 
wavers than at cruise. The main reasons limiting fur- 
ther inroads of CFD into the design process are: (1) 
the reliability of turbulence 1 models and (2) the time 
and expense 1 of the numerical simulation. Because of 
the prohibitive resolution requirements of direct sim- 
ulations at high Reynolds numbers, transition and 
turbulence modeling is expected to remain an issue 
for the* near term. 1 The focus of this paper addresses 
the latter problem by attempting to attain optimal 
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efficiencies in solving the governing equations. Typi- 
cally current CFD codes based on the use of multigrid 
acceleration techniques and multistage Runge-Kutta 
time-stepping schemes are able to converge lift and 
drag values for cruise configurations within approxi- 
mately 1000 residual evaluations. More complexity in 
the geometry or physics generally requires many more 
residual evaluations to converge, and sometimes con- 
vergence cannot be attained. An optimally convergent 
method is defined 2 ' 5 as having textbook multigrid effi- 
ciency (TME). meaning the solutions to the governing 
system of equations are attained in a computational 
work which is a small (less than 10) multiple of the 
operation count in the discretized system of equations 
(residual ('valuations). Thus, there is a potential gain 
of more than two orders of magnitude in operation 
count reduction if TME could bo attained. 

In this paper, a distributed relaxation approach to 
achieving TME for Reynolds-averaged Navier-Stokes 
(RANS) ('(illations is discussed along with the founda- 
tions that form the basis of this approach. Because the 
governing equations are a set of coupled nonlinear con- 
servation equations with discontinuities (shocks, slip 
lines, etc.) and singularities (flow- or grid-induced), 
the difficulties are many. The TME methodology in- 
sists that each of the difficulties should be' isolated, 
analyzed, and solved systematically using a carefully 
constructed series of model problems. An important 
aspect of the distributed relaxation approach is a sep- 
arate' treatment of each of the factors (elliptic and 
hyperbolic) constituting the system of partial differ- 
ential equations. Another distinguishing aspect of the 
approach is that these factors are treated directly for 
steady-state flows rather than through pseudo-time 
marching methods; time-dependent, flow solvers can be 
constructed within this approach and in principle are 
simpler to develop than steady-state solvers. An ex- 
tensive list of envisioned difficulties in attaining TME 
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for CFD simulations, along with possible solutions, are 
discussed elsewhere.' 1, This paper also summarizes re- 
cent progress towards the attainment of TME in basic 
CFD simulations. 

Foundations for Textbook Multigrid 
Efficiency 

The basic framework for TME solvers is full multi- 
grid (FMG) algorithms. 2,3,0 ' 9 In FMG algorithms, 
the solution process is started on a very coarse grid 
where the computational cost of solution is negligible. 
The course-grid solution is then interpolated to the 
next fine grid to form an initial approximation. Few 
multigrid full approximation scheme (FAS) cycles, or 
possibly just one, are performed next to obtain ail im- 
proved fine-grid solution approximation. Then, the 
process proceeds to finer grids until the solution on 
the target finest grid is achieved. 

In solution of highly nonlinear problems, a good ini- 
tial guess is important. A general way to obtain such 
an initial guess is by continuation, in which the so- 
lution to the target problem is approached through 
the solutions of a sequence of parameterized problems. 
Usually the problem starting the continuation process 
is easy to solve, and difficulty gradually increases with 
control parameter approaching the target value: this 
continuation process can often be integrated into an 
FMG solver. For example, with viscosity as the control 
parameter, at the coarse grids more artificial viscosity 
can be used, then gradually be taken out as the algo- 
rithm proceeds to finer levels. Such FMG continuation 
is often natural because larger numerical viscosity is 
introduced on coarse grids, even without aiming at 
continuation. 

A version named A-FMG algorithm provides the 
device needed for optimal adaptive local refinement. 
Efficient multigrid solvers based on this approach have 
been demonstart ed. 10 

Tlu 1 objective of FMG algorithms (and TME meth- 
ods in particular) is fast convergence to the solution of 
the differential equations, not necessarily fast asymp- 
totic residual convergence. The natural solution toler- 
ance is the discretization error defined as the difference 
between the exact solutions of discrete* and differential 
problems. Thus, the quality of a solution approxima- 
tion on a given grid can be measured by the relative 
magnitude of algebraic errors in comparison with the 
discretization error level. The algebraic error is defined 
as the difference between the exact and approximate 
solutions of the discrete problem. On any grid in an 
FMG algorithm, we expect the algebraic errors after 
few multigrid cycles to he always loss than the dis- 
cretization error. 

On the other hand, a fast residual convergence is 
considered as an important- monitoring tool. In many 
practical cases, it is possible to develop a solver ex- 
hibiting fast residual convergence rates without com- 


promising TME. Note however that sometimes the 
quality of the target-grid solution can be much im- 
proved by double discretization methods applying for 
relaxation a different scheme than that used in cal- 
culating residuals transferred to the coarse grid; zero 
target-grid residuals might not be the aim in this case. 

Standard multigrid methods efficient for elliptic 
problems separate 1 the treatment of oscillatory and 
smooth error components. The 1 former are efficiently 
reduced in single-grid iterations (relaxation); the lat- 
ter are well approximated on coarse 1 grids anel, hence, 
eliminated through the 1 eoarse-grid correction. The dif- 
ficulties associated with extending TME for solution of 
the RANS equations relate to the fae*t. that these equa- 
tions arc 1 a system of coupled nonlinear equations that 
is not, even for subsonic Mach numbers, fully ellip- 
tic, but contain hyperbolic partitions. The efficiency 
of classical multigrid methods severely degrades for 
nonclliptic problems because some smooth characteris- 
tic components cannot be adequately approximated on 
coarse grids. 11 13 The characteristic components are 1 
much smoother in the characteristic directions than in 
other directions. To be 1 efficient., a multigrid solver for 
nonelliptic problems has to adequately address three 
types of errors: (1) high-frequency error components, 
(2) uniformly smooth error components, (3) character- 
istic* error components. 

If the target discretization is strongly h- elliptic 
(or sejni-h- elliptic) one can design a local (or block- 
wise) relaxation procedure efficiently reducing all high- 
frequency error components. By definition, 2,3,8, 1 1 a 
discrete scalar (not necessarily elliptic) operator L[u] 
possesses a good measure of //-ellipticity, if the abso- 
lute value of its symbol \L{6)\ = |r“' ( ^ ) L[c ?(/7, J ) ]| is 
well separated from zero for all high-frequency Fourier 
modes. Here j — (j x . ) are the grid indexes and 

0 = (0 X , O y ,f) z )A) < |0. r |, \&y\, \0 Z | < 7T are normalized 
Fourier frequencies. High-frequency Fourier modes an 1 
the modes satisfying max(|0. f |, |0*|) > f . For sys- 

tems, the measure of /^-ellipticity is defined as the 
measure' of the 1 determinant operator. 

Coarse-grid correction is usually efficient, for uni- 
formly smooth error components. An effective reduc- 
tion of characteristic error components can be achieved 
either by designing a proper relaxation scheme 1 reduc- 
ing not only high-frequency but smooth error com- 
ponents as well (which can be done in many non- 
uniformly-elliptic cases by downstream ordering of re- 
laxation steps 11,12 ) or by adjusting coarse-grid opera- 
tors for a better characteristic-component approxima- 
tion. 

Multigrid methods efficiently reducing all the three 
aforementioned types of error have been developed for 
scalar nonelliptic operators. 13 17 Similar efficiency for 
solution of the R ANS system of differential eejuations 
can be achieved by exploiting the system factorizabil- 
ity. Fact o7-i zabihty is a property of the system deter- 
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minant to be factored as a product of simple scalar 
factors representing the (elliptic and hyperbolic) par- 
titions contributing to the target differential system. 

For the subsonic compressible Euler equations, the 
first TME solvers exploiting factorizabilitv of the sys- 
tem have been developed by Ta’asan. 18 20 New canon- 
ical variables have 1 been introduced, and in these vari- 
ables. the Euler system of equations has a block upper 
triangular form with the main diagonal blocks con- 
sisting of the basic components of the system. 21 The 
main disadvantage of this formulation is that it is not 
easily generalized to viscous and unsteady problems, 
especially in three dimensions. 

Another approach toward achieving TME for so- 
lution of the Euler and incompressible Navier-Stokes 
equations 22 " 24 is based on the pressure-equation for- 
mulation which effectively separates elliptic and hy- 
perbolic factors of the system. This formulation 
has been extended to generalized coordinates. 2> 2f) 
Grid-independent convergence rates have been demon- 
strated for inviseid flow around airfoils and for viscous 
incompressible flow past a finite flat plate and a non- 
lifting airfoil. This approach has met difficulties in 
generalizing to viscous compressible flows. 

A more general approach is the distributed relaxation 
method. 2 ’ 3 ' 12,27 29 The details of the method are pre- 
sented in the next section. An important feature of 
this approach is that the separation of different, fac- 
tors constituting the system determinant occurs only 
in computing iqxlat.es at the relaxation stage of an 
outer multigrid solver: the original coupled (conser- 
vative) equations are always used to compute resid- 
uals. This feature allows a lot of freedom in relax- 
ation scheme design since different schemes may be 
applied to different flow regions. The distributed re- 
laxation approach reduces the problem of relaxing a 
complicated system of discretized coupled differential 
equations to relaxation of scalar factors constituting 
the system determinant. Note that relaxation schemes 
for the scalar factors may include separate multigrid 
solvers. Usually, distributed relaxation can be ap- 
plied throughout the entire domain having the full 
effect away from discontinuities (shocks, slip lines) in 
the regular (smoothly varying) flow field. Some local 
relaxation sweeps should be applied in these special 
regions after (and perhaps also before) the distributed 
relaxation pass to reduce residuals. The general rule 
for efficient adaptive 1 relaxation is to apply additional 
relaxation sweeps wherever local residuals significantly 
exceed the average level characterizing the neighboring 
regular flow field. 

The distributed relaxation scheme design for the 
RANS system of equations can be significantly simpli- 
fied if the target discretization is also factorizable, i.e., 
the discrete system determinant can be represented as 
a product of discrete scalar factors, each of them ap- 
proximating a corresponding factor of the determinant. 


of the differential R ANS equations. In fact, since dis- 
tributed relaxation is applied only for solution updates 
in a relaxation sweep, the factorizabilitv property is 
only required for the principal linearization opera- 
tor. The principal linearization of a scalar equation 
contains the linearization terms that make* a major 
contribution to the residual per a unit change. The 
principal terms thus generally depend on the scale, or 
mesh size, of interest. For example, the discretized 
highest derivative terms are principal on grids with 
small enough mesh size. For a discretized system of 
differential equations, the principal terms are those 
that contribute to the principal terms of the system 
determinant. If the principal linearization is discretely 
factorizable and efficient relaxation schemes for the 
corresponding discrete scalar factors are available, dis- 
tributed relaxation efficiently reduces high-frequency 
and characteristic error components as well. 

For nonfact or izable (but /i-elliptic) discretizations of 
the RANS equations, the general scheme for relaxation 
updates should include two different passes: ( 1 ) Di- 
rect relaxation of the target discrete scheme that is 
efficient for high-frequency error reduction and (2) dis- 
tributed relaxation based on reasonable discretizations 
of the scalar factors of the differential system determi- 
nant. (these discretizations are not derived from the 
target discrete system) that eliminates characteristic 
error components. 

Distributed Relaxation 

The system of time-dependent compressible Navier- 
Stokes equations can be written as 

0 t Q + R(Q) = 0, (1) 

when' the conserved variables are Q = 
(////, pv, pu\ p, pE ) 7 , representing the momentum 
vector, density, and total energy per unit volume, 
and R(Q) is a spatial divergence of a vector function 
representing convection and viscous and heat transfer 
effects. In general, the simplest form of the differential 
equations corresponds to nonconsorvative equations 
expressed in primitive variable's, here taken as the set 
composed of velocity, pressure, and internal energy, 
q = ( u . r. ic. p, e ) 7 . For a perfect gas, the primitive 
and conservative variables are connected through the 
following relations 

P = ( 7-1 )/*7 

e = E - i (^u 2 + v 2 + w 2 ^j , 
r 2 = 7 p/p, 

where c is the speed of sound and 7 is the ratio of 
specific heats. 

The time-dependent nonconservative equations are 
found readily by transforming the time-dependent con- 
servative 1 equations. 
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7)q [Ot Q + R] = 0, 

^ < 1 + Iq R ' = 

where is the Jacobian matrix of the transforma- 
tion. For steady-state equations, the time derivative 
is dropped. In an iterative procedure, the correction 
= q' 14 ’ 1 — q". where u is an iteration counter, (*an 
he computed from the equation 




( 2 ) 


where L is a linear operator containing the viscous and 
inviseid terms of the nonconservative operator that 
are principal at the scale h. Thus, a good correction 
is expected away from discontinuities in the regular 
(smoothly varying) flow field. 

Usually, the principle linearization operator L is de- 
rived under the assumption of solution smoothness 
that requires the magnitudes of all solution differences 
to be smaller than the solution magnitude. The ap- 
proximation inherent in this principal linearization can 
be illustrated by using the nonlinear convection ('([na- 
tion. 


the entire system associated with L. In relaxing scalar 
factors, the changes introduced in the "ghosU vari- 
ables Sw (the variables Sw an* “ghost** because* they 
need not be* explicitly used in computations) during 
relaxation are distributed, with the pattern of distri- 
bution matrix M, to tin* primitive variables. To obtain 
the optimal (textbook) efficiency, relaxation of each 
factor should incorporate* the essential part of an ef- 
ficient multigrid solver for its corresponding operator: 
sometimes this essent ial part is just the relaxation part 
of that solver, sometimes this may even be an entire 
separate multigrid solver applied at souk* proper sub- 
domains. 

Incompressible Navier-Stokes Equations 

The steady-state incompressible Navier-Stokes 
equations can be written as 

Qv u + S7P = 0. 

V ' u = b, 

where u = (u. c. tr) 1 is the velocity vector and Q„ = 
u ■ v ~ ^ i s a convection-diffusion operator. (Q = Qu 
denotes the particular case with zero {u = 0) physical 
diffusion.) Under the solution smoothness assumption 
the principal linearization operator is given by 


N(u) — udj.u = /. 

A full linearization for a correction Su results in 

‘——Su = Sudj.u + udfSu — f — N(u). 

TIk* principal linearization of this correction equation 
at scale h is 


Q v 0 0 0* 

0 Qv o o f , 

0 0 Qv 0 Z 

(% Of, d : 0 


(3) 


det L = ~~Qi -X (4) 


udj.Su = f — N (a ) , 

where the term S tidbit can be neglected as h — > 0 as- 
suming that d } ‘ u is bounded. This approximation is 
also termed Picard iteration, which is exact for the lin- 
ear case. Note that on coarser grids, the term Sud f *u 
may not be so small. The FMG algorithm plays a very 
important role in preventing fine-grid initial approxi- 
mations with large high-frequency algebraic errors vi- 
olating the smoothness assumption. 

While significantly simplified by retaining only prin- 
cipal terms, the system (2) is still a set of coupled 
equations containing elliptic and hyperbolic compo- 
nents. Therefore, collective Gauss- Seidel relaxation of 
L is not often effective, and factorizability of L must 
be exploited. The distributed relaxation method re- 
places riq in (2) by Mriw, so that the resulting matrix 
L M becomes lower triangular. The diagonal elements 
of LM art* composed ideally of the separable factors 
of the matrix L determinant. These factors are scalar 
differential operators of first or second order, so their 
efficient relaxation is a much simpler task than relaxing 


where the coefficients ( u, e, w) in Q u are computed 
from the previous solution approximation and fixed 
during each distributed relaxation step. An appropri- 
ate matrix M is 


M = 


1 0 0 -0 T 
0 1 0 -d„ 

0 0 1 -<j. 

0 0 0 Q, 


yielding the lower triangular operator 


LM = 


Q v 0 0 0 

0 Qr 0 0 

0 0 Q„ 0 

0, d„ 0 - -A 


(5) 


( 0 ) 


Euler Equations 

The conservation form for the Euler equations is 
given bv (1) with 


R(Q) = 0,F(Q) + 0,G(Q) + fl : H(Q), (7) 
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F(Q) = 


( f»r + V \ 

f)UV 

paw 
pu 

\ puE + up ) 


G(Q) = 


/ puv \ 

f)V 2 + P 

pvw 

pv 

V pt’E + vp ) 


/ 


H(Q) 


\ 


paw 
pvw 
pw 1 + p 
pw 

y pwE -f wp } 

The Jacobian matrix of the coiiservative-to- 
nonconservative transformation is defined as 


dq 

OQ 


( i Ip o 

o l Ip 
0 0 
0 0 
\ 0 0 

/ 1 0 0 

(j 1 0 

0 0 1 

-u —v -w 

V — w -v —w 


0 

0 

i Ip 
o 
o 


0 

0 

0 

0 


1 


0 \ 

0 
0 
0 

1 Ip j 

0 \ 

0 

0 

1 

f + i J 


— u 

— V 

—w 
u’ + r- + 


In the l egions of smoothly varying solution, the prin- 
cipal linearization of the noncouservative operator is 


Q 0 0 

0 Q 0 

0 0 Q 

pc' 2 0 r p( :2 0„ pc 2 d- 

L T«- T 3 » 


0 

0 


\ 0 ' 

l -d- 0 
/' - 

Q o 

0 Q 


( 8 ) 


The determinant of the matrix operator L is 
det L = Q 3 [Q- - c 2 A] , 


(9) 


where A is the Laplace operator, and Q 2 - < J A repre 
sents the full-potential operator. 

A possible distribution matrix M is given by 


M = 


1 0 0 - L dx 


0 1 0 


Du 


0 0 1 ~jd z 

0 0 0 Q 
0 0 0 0 


0 

0 

0 

0 

1 


and 


LM = 


Q 0 0 0 

0 Q 0 o 

0 0 Q 0 

pc 2 dj. P< :1 d u pr 2 d z Q 2 - r* 2 A 

( ^0.r -0, 


- — A 


0 

0 

0 

0 

Q 


Compressible Navier-Stokes Equations 

The conservative compressible Navier-Stokes equa- 
tions are formulated in the form (1) with R(Q) defined 
in (7). 


F(Q) - 


G(Q) = 


( pu 1 + p — 2pO x u - A( V ' u ) \ 
par - p{Oyv + Oyti) 
paw — p(0, v w 4- 0 : a) 
pu 

y pu E 4- up — A i/.(V ' u ) — P T i ~ Kd x f J 


( pu v p { Oj. t ' 4- 0 f/ u ) \ 

pv 2 4- p - 2pO f ,v - A( V * u ) 
pV W — p (Of, W 4- 0 Z I*) 
pv 

y pvE + vp - A r(v • u) - pr -2 - h’d u f / 


H(Q) - 


puw - p(dyw 4- 0 z u) 
pvw — p(Oy W 4- 0 z v) 

pw 2 -bp - 2 pO z w - A(v * u ) 

pw 

y pwE + wp - Atr(V • u) “ l ,T -< ~ K d- ( J 


where 

n = 2udj.u 4- v(0j. v + OyU) + w(O r w 4- 0 z u ), 
r> = 2 vOyV + u(dyv 4 0 u u) 4- w(OyW 4- 9 : r). 

= 2 wO z W + u{0 r w 4- 0 z u) -b v(OyW -b 0 Z V ) . 

p and A are viscosity coefficients, and k is the coeffi- 
cient of heat conductivity. 

The corresponding noncouservative' formulation is 
given by 


( 10 ) 


.(ID 


f(u • V) - p*. - p&rs)‘ u - p(P’!/ V + 0 xz w) + jdj-p = 0. 

((u ■ v) - ^ - Jpvy) v - jS d ry u + ^ = (, ‘ 

^(u • V) ~ p A - w + 0,, : V) + yp-j) = o. 

pc 2 ( v • u) + (u • y)p + (', - 1 )(-kAc + *!>)= 0, 

^(V • u) + ( u • V)P ~ fA ( + P® = l) - 

wh(*re 

$ = //^2(9j.w) 2 + 2 (Oyi') 2 -b 2 (0 z w)~ 

+ {Q x v 4- OyU)- 4 (Oyw 4- 0 z u) 2 4- (OyW -b 0 z v)“ 
+A {Oyu -b OyV 4- 0 z w) 2 . 

Assuming solution smoothness and constant viscos- 
ity and heat conduction coefficients, the principal lin- 
earization operator L. keeping the terms principal on 
both the viscous and inviscid scales, is given by 
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0 


L - 


d - 

,■ --0 r „ 

-fa- 

~l 

y (> 

f) 


f> ’ 

~fa» 

(h - fa 


fa 

-fa 

'-fa 

- fa 

±t>- 

p 

f)C~0j 

/'CO,, 

fCO: 

Q 

^ 0. r 

fa 

C-0: 

0 


0 
0 

(1 - *, )k4 

Q* 


( 12 ) 


dot L = Q'i 


+ Q(-, :1 A + A- ) 

_g-.>K±A±ii a +(/'], 


(13) 


whore, nondimensionalizing by density and sound 
speed and applying Stokes hypothesis for the bulk 
viscosity term, the coefficients become p/p = 
/ (p Re), k = M x 7/(Re Pr), and A = A + // = /// 3, 
A/, x is the free stream Mach number, and Re and Pr 
an' Reynolds and Pranritl numbers respectively. In- 
stead of devising a suitable relaxation method for the 
complicated scalar factor in the brackets of (13). one 
can opt to a distributed relaxation partially decoupling 
the linear system associated with operator L (12). In 
particular, the distribution matrix 





l 

0 

» -fa 

0 ■ 





0 

1 

« ~fa 

0 


M 

= 


0 

0 

1 ~fa 

0 





A ch 

A 0 U 

A(T Q 

0 





0 

0 

0 l) 

1 _ 


results in 










l) 

0 

0 

0 

~ 


0 


Qh 

0 

0 

0 


LM = 

0 


(/ 


0 

(} 



P0 r 

POy 

P0:_ 

QQx±rc* A (i 

- 7)kA 


far 


T* 

-^-A 

Q : 

J 


(Id) 


(15) 


where V = p< :1 + A Q. The last two equations re- 
main coupled, requiring a block 2-by-2 matrix solution. 
This distributed relaxation scheme is still much less 
expensive than direct relaxation of matrix L requiring 
solution for a block 5-by-5 matrix. 


Relaxation of Scalar Factors 

Efficiency of the distributed relaxation schemes out- 
lined in the previous section is determined by the ef- 
ficiency of the relaxation (solution) schemes for scalar 
factors appearing at the main diagonal of the matrices 
LM. 

For uniformly elliptic operators such as Lapla- 
cian. diffusion-dominated convection-diffusion opera- 
tor, and subsonic full-potential operator many effi- 
cient relaxation techniques are available (see text- 
books ,u> ’ s,y ). For such operators, an important re- 
laxation requirement, is efficient reduction of high- 
frequency errors. All the smooth components are well 


approximated on coarse grids built by standard (full) 
coarsening: therefore, the eoarse-grid correction is ef- 
ficient in reduction of smooth errors. 

For nonelliptic and weakly elliptic factors, e.g., con- 
vection, convection-dominated convection-diffusion, 
transonic and supersonic full-potential operators, 
(smooth) characteristic components cannot he approx- 
imated with standard multigrid methods. 11 

Several approaches aimed at curing the 
characteristic-component, problem have been studied 
in the literature. These approaches fall into two 
categories: (1) development of a suitable relaxation 
scheme (single-grid method) to eliminate not only 
high-frequency error components but the characteris- 
tic error components as well: (2) devising an adjusted 
eoarse-grid operator to approximate' well the fine-grid 
characteristic error components. 

Single-Grid Methods 

D o wnstream m arching 

For hyperbolic problems, the simplest first-category 
method is downstream marching. If the c orrespond- 
ing discretization is a stable upwind discretization and 
the characteristic field does not recirculate, then down- 
stream marching is a very efficient, solver that yields 
an accurate 1 solution to a nonlinear hyperbolic equa- 
tion in just a few sweeps (a single 1 downstream sweep 
provides the exact solution to a linearized problem). 
Tilt 1 downstream marching technique 1 was successfully 
applied in solving many CFD problems associated with 
non-recirculating flows (see. e.g., 12, - :! ' a5 -- 7 ’ 28 ) i How- 
ever. if a discretization operator is not fully upwind 
(e.g.. is only upwind biased), straightforward down- 
stream marching is unstable 1 . For the schemes that 
cannot be directly marched, there are two possible al- 
ternatives (also of marching type): defect-correction 
and predictor-corrector methods. 

Defect Correction 

Let us consider a defect correction method for a dis- 
cretized hyperbolic equation 

L h n iuh = f h . h , (1C) 

with specified inflow boundary conditions 

Let Ui x j 2 be the current solution approximation. 
Tin'll the improved approximation f/ f - lW2 is calculated 
by defect-correction scheme in the following two steps: 

L The correction e,, ., 2 is calculated by solving oper- 
ator Lj with a right- hand side represented by the 
residual of (1G) computed for the current approxi- 
mation . The inflow boundary conditions for 
v are initialized with the zero value's. 

LrfViiJj ~ fi\*h ~ ^ (10 

2. The current approximation is corrected as 

« m . i 2 = u ,,.; 2 + (18) 
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The operator V' (l is called the driver operator. It is 
chosen to be easily solvable and usually less act u- 
rate than the target operator L h ; the latter can be 
very general. If the iteration converges, steps (17) 
and (18) can be repeated until the desired accuracy 
is reached. Usually the efficiency of defect-correction 
methods is quite satisfactory." ‘ -S, il even though 
in principle the convergence rate of a defect -correct ion 
method for nonelliptic oj)erat.or is normally mesh-size 
dependent.. 31 36 as explained below. 

In several papers (e.g. 32 * 37 ), authors studying the 
defect-correction method for nonelliptic problems ob- 
served a slow convergence or even a divergence 1 in some 
common error norms for the initial iterations and good 
asymptotic convergence rates afterward. This behav- 
ior is different from that observed in solving elliptic 
problems by the defect-correction method, where the 
asymptotic convergence rate is the slowest one. This 
nonelliptic feature is explained by some properties as- 
sociated with the cross-characteristic interaction (e.g.. 
dissipation and/or dispersion) in the operators in- 
volved in the defect -correction iterations. Specifically, 
this cross-characteristic interaction defines the pene- 
tration distanee (also termed “survival distance 1 ) 
of a characteristic component.. The penetration dis- 
tance is the distance from the inflow boundary within 
which the discrete solution of the homogeneous prob- 
lem reasonably approximates the continuous one (i.e., 
the discretization error is substantially smaller than 
the solution). 

The penetration distance of a characteristic com- 
ponent. is roughly proportional to u; ii , where 

q is the highest order of differentiation in the hy- 
perbolic operator under considerations, p is the 
discrete-operator approximation order, u; is the cross- 
characteristic frequency of the characteristic compo- 
nent, and h is the mesh size. The ratio of penetration 
distances of the operators and Tj} is an important 
factor for determining the number of defect-correction 
sweeps required to reduce the algebraic error to the 
discretization-error level or to reach the asymptotic 
convergence regime. 

When the operators L h and Lj} have the same 
approximation order (p = r). efficiency of the 

defect-correction method is optimal and mesh-size h- 
independent If however the operators L }l and L [ j have 
different, approximation orders {p and r, respectively, 
p > r), then efficiency of the defect-correction method 
is //-dependent; i.e., the maximal number of sweeps 
which might be required to reduce the algebraic er- 
ror to the discretization-error level (or to reach the 
asymptotic convergence rates) is larger on fine grids 
than on coarse grids. This is because one has to iterate 
L[\ as many times as needed to attain accuracy up to 
the L h penetration distance. The worst (largest ) ratio 
of penetration distances is obtained for characteristic 
components for which the penetration distance in the 


target p-order accurate discretization approaches the 
depth R of the computational domain in the charac- 
teristic* direction. It follows that the required numbei 


of iterations is proportional to f j 


Predictor- Corrector 

One potentially efficient but yet unexploited 
method to overcome grid-dependent convergence 1 expe- 
rienced in defect-correction iterations is the predictor- 
corrector technique. A detailed look into the defect- 
correction iteration reveals that the computational 
work distribution is unbalanced: (1) Driver operator 
iterations at locations beyond the penetration distance 
do not improve 1 the solution approximation. (2) In suc- 
cessive iterations, the solution approximation near the 
inflow boundary becomes much more accurate than in 
the interior; the computational efforts spent in this re- 
gions could be more profitably invested at the accuracy 
frontier. 

The predictor-corrector method has been exten- 
sively used for ordinary and time-dependent differen- 
tial equations; 38 * 39 however, applications for steady- 
state* nonelliptic problems have been very limited. In 
predictor-corrector schemes, the final update of the so- 
lution at a particular point is computed from the local 
solution of the target operator. The solution values 
at downstream points included in the target-operator 
stencil are predicted from the solution of the dii\ei 
(predictor) operator. In order to define? a family of 
predictor- corrector schemes, one can divide the com- 
putational domain into several time-like layers: the 
first layer contains all the grid points adjacent to the 
inflow boundary. Each next layer is composed of the 
grid points that contribute to the stencils of target op- 
erators defined at the points of the previous layer and 
do not belong to any of the previous layers. 

Now, a family of predictor-corrector schemes for 
solving the correction equation 

L ,l v h j. i ^R' l \j . 2 = f L h u n . h , (19) 

where L h is the target operator, u h j 2 is the current 
solution approximation with residual , and Vj , 

is the desired correction function, can be defined as 

PC {) : The solution of (19) is approximated by solving 

L\ }v h . h =Rl h - ( 20 ) 

This scheme 1 is identical with the defect-correction 
scheme. 


PC k : 


Recursive definition of the derived predictor- 
corrector schemes (recursion with respect, to k) 
can be done as follows: Assume the (j — l)-th lavei 
have already been finally updated in the current 
sweep. Then, to calculate new* values at the next 
j-tli layer one has to perform the following three 1 
steps: 


of r 


American Institute of Aeronautics and Astronautics Paper 2001 2570 



1 To predict values at the j - th layer with 
PCk- i scheme. 

2. To predict values at the {j + 1 )-t h layer with 
PCb _ i seh(*me. 

3. To update values at the j-t h layer by directly 
relaxing (19). 

The schemes for k = 0,1.2 have been tested for a 
linearized supersonic full-potential operator. 17 

Multigrid methods 

ReciiTulation 

Downstream marching methods an' not viable for 
problems with closed characteristics. Alternative dis- 
cretization and solution techniques should be consid- 
ered. The discretization issue becomes especially im- 
portant for flows with streamlines that do not start 
and end at boundaries, but constitute closed curves. In 
such cases, even a very small viscosity plays an impor- 
tant role in determining the flow throughout the do- 
main. The solution in the limit of vanishing viscosity 
depends very strongly on how the viscosity coefficients 
tend to zero. The propagation of information from the 
boundary into the domain is governed by the viscous 
terms no matter how small they may be. It has been 
shown 10 for both the scalar convection-diffusion prob- 
lem and the incompressible Navier-Stokes equations 
that varying cross-stream numerical viscosity (caused 
usually by varying angles between the stream and the 
grid lines; e.g., in standard upwind and upwind bi- 
ased scheme's) may prevent convergence to a physically 
realizable solution. In the most general case. it. can 
be shown that even isotropic viscosity is not suffi- 
cient for convergence to a physical solution, and one 
must actually specify a uniform viscosity. However, 
for the homogeneous problems there air several indi- 
cations 40 * 41 (though no proof) that isotropy suffices. 

To obtain a discretization scheme that exhibits the 
appropriate physical-like behavior for vanishing viscos- 
ity, one must either add sufficient explicit isotropic* 
viscosity that will dominate the anisotropic numeri- 
cal viscosity of the convection operator, or else 1 derive 
a discrete convection operator with numerical viscos- 
ity satisfying the condition of isotropy. An upwind 
isotropic- viscosity discretization has been derived. 11 

One general approach to the algebraic solution 
of nonelliptic equations with closed characteristics is 
to apply a multigrid method with an overweighted 
upwind- biased residual restriction. Efficient multigrid 
solvers for recirculating convection equation have al- 
ready been demonstrated. 41 42 This approach is well 
combined with the distributed relaxation method for 
the RANS ('(illations, because* within a distributed 
relaxation sweep a multigrid solver with optimal over- 
weighting can be applied to a separate scalar nonellip- 
tic equation with closed characteristics. 


Another solution approach is to apply some gen- 
eral techniques to approximate indirectly smooth char- 
acteristic components. Among helpful techniques 
are recombination of iterants. cycles with high in- 
dexes, and implicit alternative-direction relaxation. 
Recombination of iterants (solution approximations 
on different stages of a multigrid algorithm) at each 
grid level eliminates several (number of iterants mi- 
nus one) error components, those, more* specifically, 
that are most prominent in the residual function. 
Making increasingly many course-grid iterations per 
each fine-grid iteration, cycles with high indexes 
solve the characteristic-component problem on coarser 
grids. Implicit alternative-direction relaxation simu- 
lates downstream marching in the regions with open 
characteristics and efficiently transfers information in 
the regions with characteristics closely aligned with 
the grid. Theoretically, each of these methods can- 
not completely resolve the problem of poor coarse-grid 
correction to tin* fine-grid smooth characteristic error 
components. The problem already manifests itself in 
two-level algorithms with any type of local relaxation. 
On fine grids the number of problematic error compo- 
nents may increase, and many cycles may be needed to 
collect, the necessary number of the fine-grid iterants to 
exclude all the troubling error components. However, 
it has been shown experiment ally 4 ** that a combination 
of implicit alternative-direction defect-correction type 
relaxation, recombination of iterants on all the levels, 
and W-cycles can result in a relatively efficient multi- 
grid solver for recirculating flow problems on practical 
grids. 

Full- Po t c n tial ( )perato r 

The full-potential operator is a variable type oper- 
ator. and its solution requires different procedures in 
subsonic, transonic*, and supersonic regions. In deep 
subsonic regions, the full-potential operator is uni- 
formly elliptic and therefore standard multigrid meth- 
ods yield optimal efficiency. When the Mach num- 
ber approaches unity, the operator becomes increas- 
ingly anisotropic and, because smooth characteristic 
error components cannot be approximated adequately 
on coarse grids, classical multigrid methods severely 
degrade. In the deep supersonic regions, the full- 
potential operator is uniformly hyperbolic with the 
stream direction serving as the time-like direction. In 
this region, an efficient solver can be obtained with a 
downstream marching method. However, downstream 
marching becomes problematic for the Mach number 
dropping towards unity, because the Courant number 
associated with this method becomes large. Thus, a 
special procedure* is required to provide ail efficient 
solution for transonic regions. A possible local pro- 
cedure 1,1, ir ** 1 l,t0 is based on piecewise semicoarsening 
and some rules for adding dissipation at the* coarse- 
grid levels. 
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Fig. 1 Staggered arrangement of primitive vari- 
ables for Navier-Stokes discretization. 

A similar technique can be used to construct ail 
efficient marching-free multigrid solver for convection- 
dominated equations. This method 16 employs a col- 
ored relaxation scheme and is very attractive for mas- 
sive parallel computing. A highly parallel multigrid 
solver for the supersonic full-potential operator may 
be obtained bv methods similar to the wave/ray multi- 
grid. 1 1 

Discrete Equations 

Traditional Factorizable Schemes 

As mentioned in the introduction, factorizability of 
a target discrete scheme significantly simplifies the 
distributed relaxation design. The main obstacle in 
this case to efficient solution is that the discretiza- 
tions obtained for the scalar factors in the discrete 
determinant are not always convenient. The search for 
new factorizable discretization schemes (see Summary 
of Recent Progress below) is chiefly motivated by the 
need to derive' discrete schemes with the resulting dis- 
cretizations of scalar factors satisfying some desired 
properties (e.g., correct alignment with the physical 
anisotropies, compactness, availability of an efficient 
relaxation scheme*, etc.). 

Staqgervd-Grid Discretization for N avicv-Stokc.s 
Equations. 

The staggered-grid discretization dating back to the 
mid GO’s 4 1 46 is one of the first factorizable discretiza- 
tions for incompressible flow equations. Compress- 
ible flow discretizations with a staggered arrangement 
of variables have also been studied 3 * 2 '* 17 A usual 
placement of primitive variables in two dimensions is 
depicted in Figure 1. With this staggering, (a) the 
off-diagonal first derivatives in (3), (8), and (12) can 
be approximated as short central differences; (b) the 
second derivatives in (12) can be compositions of cor- 
responding central first derivatives; (c) the convection- 
diffusion operators. Q „ , (“an be approximated by am 
proper discretizations Qj. For discrete factorizability, 

0 of 


it is important to have the same discretization for each 
of the (^-operators in the momentum equations; the 
convection-diffusion operators in other equations can 
be different. The convection terms in the momentum 
and energy equations are usually upwind or upwind- 
biased; for simplicity, below we assume that all these 
terms are the same. With such differencing, the dis- 
crete schemes mimic the factorizability property of the 
differential equations, and the discrete* system deter- 
minants can be factored as det L h = (Q$) 2 ± h (incom- 
pressible Xavier-Stokes) or det L * 1 — Q ~ 

c~A^) (compressible Euler), where? in three* dimen- 
sions is the* seven-point /i-Laplacian, Q h is an upwind 
or upwind biased discretization of the convection op- 
erators in the momentum and energy equations, Q 1 
is a convection-operator discretization for the pressuie 
term in the fourth equation of (8), hence Q h Q h -c 2 
is a discrete approximation to the* full-potential op- 
erator. The discrete determinant computed for the 
compressible Navier-Stokes equations is similar to the 
differential determinant (13). 

The discrete distribution matrices follow directly 
from the continuous matrices (5). (10), and (14). The* 
short central differences are used for the approxima- 
tion of all the off-diagonal first derivatives; the con- 
vection parts in the ^-operators an* the same as those 
in the momentum equations. The resulting products 
L h M h an* similar to those for the continuous prob- 
lems with the* main diagonals composed of the factors 
of the discrete determinants. 

Distributed-relaxation solvers have been success- 
fully applied to the staggered-grid discretization 
schemes for subsonic compressible 2 ' and incompress- 
ible 12,28 flow problems. 

In computing the Euler system of equations, the 
main disadvantages of the staggered-grid scheme relate 
to the discrete stencil of the full-potential operator. 
For subsonic* flow problems, the downwind differencing 
applied for the* Q h term results in a full-potential- 
operator stencil that is somewhat wide* (because of 
the Q h Q h term) and poorly aligns with tlu* physical 
(cross-stream) anisotropies in approaching the tran- 
sonic regime. For supersonic* flow, where the problem 
is purely hyperbolic, the stencil is not fully upwind 
(even if the Q h term is upwind differenced) implying 
more involved marching schemes. 

Recently, a new approach to building discretization 
schemes that allows any desired differencing for the 
full-potential factor of the system determinant, with- 
out compromising the scheme factorizability has been 
discovered. This approach is discussed subsequently 
in application to central collocated-grid discretizations 
(see also 18 ), but it applies to staggered grids as well. 

15 
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Collocated- Giid Discretizations for the Eider 
Equations. 

Another example of a factorizable scheme is a 
eollocated-grid scheme with the second-order central 
differencing for the off-diagonal first derivatives in (3), 
(8). and (12). The convection operators in the mo- 
mentum and energy equations are again upwind or 
upwind biased: the differencing of the convection term 
applied to the pressure may alternate from down- 
wind (downwind-biased) in subsonic mode to upwind 
(upwind-biased) in supersonic mode. 

A typical difficulty associated with this type of 
schemes is a poor measure of //-ellipticity in the dis- 
crete approximation for the full-potential factor of the 
system determinant. To be more specific, let us define 
the eollocated-grid discretization L h of the matrix op- 
erator (8) as 


0 1 

& 0 
0 . 
V 0 

0 Q h 

( 21 ) 

where the discrete derivatives, d'j h , d 2h , dV l . in all off- 
diagonal terms an* the wide (with mesh spacing 2//) 
second-order-accurate central- differencing approxima- 
tions. All the diagonal terms, Q h , except Q h in the 
fourth equation, are discretized with the same* second- 
order-accurate upwind (or upwind- biased) discretiza- 
tion scheme. In the subsonic regime (|u|~ = Fr + v 2 -f 
tc 2 < c 2 ) . the Q^-tenn is discretized with a second- 
order-accurate downwind (or downwind-biased) dis- 
cretization. The determinant, of the matrix operator 
L h is given by 


Q h 

o 

0 

'—0~ h 


0 

Q" 

o 


0 

0 

Q" 

pc-ff: 1 ' 

dpi> 


((?')'' [Q h Q h - c 2 A 2, ‘] . (22) 


when* A J/| is a wide discretization of the Laplace 
operator. The full-potential-operator approximation 
appearing in the brackets has two major drawbacks: 
(1) For slow velocities (|u| <$: c), the discrete oper- 
ator is dominated by the non-fr-elliptic wide Lapla- 
cian, and efficiency of any local relaxation severely 
degrades. (2) For near-sonic regimes (the Mach num- 
ber M ~ |u|/c % 1), the discrete operator stencil does 
not reflect the physical anisotropies of the differential 
full-potential operator; the discrete operator exhibits a 
very strong coupling in the streamwise direction, while 
the differential operator has strong coupling only the 
cross-stream directions. 

Several approaches to cun* the lack of ft-ellipticit.y 
(mainly in applications to incompressible-flow equa- 
tions) have been proposed in the literature (e.g.,' 19,5 °). 
Some of the approaches are associated with introduc- 
tion of additional terms increasing the measure of 


//-elliptieity in the system of equations, others pro- 
pose averaging (filtering) spurious oscillations. The 
problem of wrong anisotropies in the full-potential- 
operator has not been sufficiently investigated. In 
two dimensions, it is possible to construct a discretiza- 
tion that satisfies the following properties: (1) At low 
Mach numbers, the discretization is dominated by the 
standard (with mesh spacing h) //-elliptic Laplaeian. 
(2) For the transonic Mach numbers, the discretiza- 
tion tends to the optimal discretization 12, J7 :i0 for the 
sonic-flow full-potential operator. (3) For supersonic 
Mach numbers, the discretization heroines upwind 
(upwind-biased) and can be solved by marching. The 
problem of constructing a good high-order discretiza- 
tion for the transonic f nil- potential operator in three 
dimensions is still open. 

Non-Factorizable schemes 

The majority of discrete schemes in current use. es- 
pecially for compressible flow but also more recently 
for incompressible flow, art* based upon a flux-splitting 
approach. The basis of this approach is the solution of 
the Riemaini problem (i.e., the time evolution of two 
regions of flow initially separated bv a diaphragm) ap- 
plied on a dimension by dimension basis. This method- 
ology has enabled the robust treatment of flows with 
strong shocks and complex geometries. However, the 
derived schemes are not discretely factorizable, except 
in one dimension. 

These discrete equations have always been solved us- 
ing collective relaxation (or pseudo-tirne-steppingj in 
multidimensional multigrid algorithms. A better effi- 
ciency should be realizable with a relaxation scheme 
that efficiently reduces both the high-frequency and 
characteristic error components. Such a scheme should 
combine two different relaxation methods: (1) A local 
relaxation scheme 1 treating directly the conservation 
equations and reducing the* high-frequency error com- 
ponents. (2) For reduction of characteristic error com- 
ponents, a defect-correction (or predictor-corrector) 
method with a factorizable driver (predictor). This 
approach has not been tried as yet. 

Boundary Conditions and 
Discontinuities 

Boundaries and discontinuities introduce some ad- 
ditional complexity in distributed relaxation. The 
determinant of L M is usually higher order than the 
determinant, of L. Thus, as a set of new variables, 
would generally need additional boundary condi- 
tions. In relaxation, because the ghost variables can be 
added in the external part of the domain, it is usually 
possible to determine suitable boundary conditions for 
<5w that satisfy the original boundary conditions for 
the primitive variables. Examples are available 28 for 
incompressible 1 flow with entering and no-slip bound- 
aries. However, to construct such extra boundary 
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conditions may bo difficult and/or time-consuming in 
general. In addition, enforcing these boundary con- 
ditions causes the relaxation equations to be coupled 
near the boundaries, not decoupled as they art* in the 
interior of the domain. 

Thus, near boundaries and discontinuities, the gen- 
eral approach’* 5 is to relax the governing equations 
directly in terms of primitive variables. Several extra 
sweeps of robust (but possibly slowly converging) re- 
laxation. such as block-Newton-Kacmarcz relaxation, 
can be made in those special regions after (and perhaps 
also before) the distributed relaxation pass to reduce 
residuals to the average level characterizing the regu- 
lar flow field. The additional sweeps will not seriously 
affect the overfill complexity because the number of 
boundary and/or discontinuity points is usually small 
in comparison with the number of interior points. An 
example of TME using local relaxation at shocks and 
boundaries and distributed relaxation over the rest of 
the domain has been demonstrated. 29 

Analysis 

As mentioned previously, it is important in attain- 
ing optimal efficiency to understand all the difficulties 
that present themselves in application. Analysis meth- 
ods are quite helpful in this regard, and the main tools 
are discussed below. In iterative methods solving el- 
liptic problems, the main mechanism of convergence is 
damping of error components. In solution of hyper- 
bolic scalar equations, there is another very important 
convergence mechanism: the downstream evolution of 
the error components. In the presence of this addi- 
tional mechanism, the accuracy first achieved near the 
inflow boundary and is then propagated into the inte- 
rior of the domain. 

The recognition of this additional convergence mech- 
anism urges modifications in the standard analysis 
developed for elliptic problems. Basically, one can 
distinguish four types of analysis applied to nonel- 
liptic problems: (1) standard linear-algebraic matrix 
analvsis, (2) modified zero- mode-exclusion full-space 
Fourier mode analysis, (3) half-space analysis of the 
first differential approximations (FDA). 11,,,l " >2 and (4) 
discrete half-space analysis. Briefly, the first differen- 
tial approximation (also called modified equation) to 
a difference operator on a grid with mesh size h is the 
Taylor expansion of this difference operator in terms 
of h truncated to the first terms including the least 
nonzero power of h. The quality of an analysis applied 
to nonelliptic problems is determined by how well the 
analysis handles the characteristic components. 

Matrix analysis 

The* most general and precise analysis methods are* 
the* linear-algebra matrix analysis methods applied to 
the* corresponding linearized problem. This analysis 
considers the* difference operators without assump- 


tions about the solution and boundary conditions. It 
can be applied t.e> variable-coefficient problems as well. 
This analysis was found very useful for analyzing one- 
dimensional problems. However, the enormous compu- 
tational complexity of this analysis makes it not viable* 
fen* multidimensional problems. Although, the analysis 
complexity can be reduced considerably by assuming 
Femrier mexles in two of the three spatial directions. 

Modified full-space Fourier mode analysis 

The modified full-space Fourier mode analysis is a 
modification of the* standard full-space Fourier mode* 
analysis excluding from the* consideration all the* zero 
modes (the models with vanishing symbols). 9 It is the 
simplest and most popular type of analysis (e.g.. see* 
applications in 22,33 ). This analysis estimates only the 
amplification (damping) factor. Its inherent disadvan- 
tage is the inability to take the influence of the inflow 
boundary into account.. This explains its failure in 
describing the downstream error evolution. However, 
the modified full-space analysis can also be useful for 
analyzing the effect of forcing terms. 

FDA half-space analysis 

The* FDA half-space analysis is a relatively simple 
and efficient tool for analyzing the effect of the inflow 
bounelarv. Examples e)f applications of this analysis 
an* available. 1 1 * 12,51,52 The first differential approxi- 
mations are* ce>nsidere*d on a half-space* including cut 
by an inflow boundary. The boundary conditions are 
represented by one Fourier moele at a time*. The* F DA 
analysis provieles a gooel qualitative description of the 
downstream error evolution. This analysis focuses 
on characteristic components and. therefore, consid- 
ers homogeneous problems. Note that a combination 
of the FDA analysis with the modified full-space anal- 
ysis e-an provide* a good insight for nonhoinogeneous 
problems as well. The disadvantages of this analv- 
sis are the inability t.e> provide quantitative estimates, 
to analyze the effect of different boundary condition 
discretizations, and to address the asymptotic conver- 
gence rate. 

Discrete half-space analysis 

The* discrete? half-space analysis 20, ,i,I ' 3t> considers the 
discretizations in their exact form rather than their 
differential approximation, while the boundary data 
are represented by a Fourier component. This analy- 
sis translates the original multidimensional problem 
into a ono-dimensie>nal discrete problem, where the 
frequencies of the boundary Fourier components are 
considered as parameters. To regularize the half-space 
problem, the solution is not allowed to grow faster than 
a polynomial function. This tool is very accurate: it 
can be used to explain in detail many phenomena ob- 
served in solving nonelliptic equations and prov ides a 
close prediction of the actual solution behavior. 

The one-dimensional solution obtained in the dis- 
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cre*t.e half-space analysis has two different representa- 
tion forms: (1) away from the boundary, the solution 
is defined as a linear combination of a finite number of 
analytical components; this region is called the analyt- 
ical representation reyion : (2) in the region adjacent to 
the inflow boundary, the solution is defined pointwise; 
this zone is referred to as the pointwise representa- 
tion, reyion. With each further iteration described 
by the analysis, the pointwise representation region 
penetrates by a finite number of mesh sizes into the 
interior. By using these representations, the compu- 
tational complexity of the analysis becomes much less 
than that associated with the one-dimensional matrix 
analysis. In the asymptotic regime, when the point- 
wise representation zone covers all the domain, this 
analysis becomes a discrete one-dimensional matrix 
analysis of the multidimensional problem. 

The discrete' half-space' analysis provides a quantita- 
tive description of the approximate solution: it predicts 
the convergence rate for each iteration and the asymp- 
totic convergence rate. It can be easily adjusted to 
analyze the global effect of any local discretization of 
the inflow boundary conditions. This adjustment can 
be done just by widening the initial pointwise repre- 
sentation region at the* inflow boundary. If necessary, 
the analysis can take into account the influence of 
the discretized outflow boundary conditions as well. 
Generally, this discrete half-space analysis treats com- 
pletely both mechanisms of convergence, damping and 
downstream evolution of errors, associated with nonel- 
liptic problem solvers. 

Summary of Recent Progress 

Tilt' first TME solver applying the distributed re- 
laxation approach for solution of a free-stream in- 
compressible Navier-Stokes equations has been demon- 
strated long ago. 12 Recently, TME has been achieved 
for high-Reynolds-number incompressible wake flow 
and the boundary layer How associated with a fi- 
nite* Hat plate. 28 An initial extension of this work 
to compressible (subsonic) viscous How has also been 
completed. 21 In all these calculations, a staggered ar- 
rangement of variable's on Cartesian grids has been 
used. With distributed relaxation, the system of equa- 
tions has been decomposed (i.e., factored) everywhere*, 
except, near boundaries where the equations remained 
coupled. Two-dimensional FMG solvers with just one 
multigriei cycle 1 per griel Jewel and a total computa- 
tional work equivalent, to about 10 target-grid residual 
('valuations converged the drag to the discretization 
accuracy. 

Recently, a new mult idimensional factorizable 
scheme for the Euler equations has beeni developed™ 
for Cartesian coordinates and extended through gener- 
alized coordinates to external lifting Hows around air- 
foils with both suberitical and supercritical freest ream 
Mach numbers. 04 -'™ This scheme is the first fiux- 


difference-splitting discretization factorizable in mul- 
tiple dimensions. The starting point for the scheme 
is the* first-order discretization of the flux-difference 
splitting scheme of Roe. Correction terms arc' added in 
the' form of mixed derivatives to make the scheme both 
second-order accurate' and discretely factorizable. The' 
resulting scheme is second-order accurate and compact 
in comparison to other scheme*. Discrete 1 fartorizabil- 
ity is achieved by using some' non-standard wide ap- 
proximations for spatial derivative's to ensure t hat the 
identities 

,r y — 0 X y0 X y , 

r dy — $jr if 1 

dy ij Of O y IJ 0 If 

are 1 satisfied on the discrete level. The 1 determinant 
of the 1 resulting schemes is compost'd of an upwind 
differe'iiced convection factor and an //-elliptic approx- 
imation for the full-potential factor. The 1 distributed 
relaxation is possible 1 by using a loft and right distri- 
bution matrix, although this has not been applied as 
yet. 

In numerical tests performed for this scheme, the 1 
multigrid solver employed alternating-direction col- 
lective Gauss- Seidel relaxation. The alternating- 
direction relaxation is necessary since the* full-potential 
factor is not. separately treated. Computations for sub- 
sonic and transonic channel flows with essentially grid- 
independent convergence rates have been presented. 54 
Grid- independent convergence rates have also been 
attained for a flow with stagnation points. 55 The 
subsonic-flow convergence rate's observed in multigrid 
V- cycles were quite fast (about 0.3 per cycle) and only 
slightly grid dependent. The rates somewhat deteri- 
orate in transonic/superscmir computations. Further 
developments of this scheme are* presented in two pa- 
pers at this conference. 55, 50 The scheme' applies at low 
Mach numbers although it has yet to be extended to 
viscous Hows. 

Another approach to building factorizable schemes 
with suitable discretizations for scalar factors has 
been explored in papers of the second and third au- 
thors. 48 * The* approach is based on a rollocated-grid 
scheme with a mechanism that allows one to improve 
the /i-ellipticity measure by obtaining any desired dis- 
cretizations for t,h<* full-potential factor of the* system 
determinant without compromising the discrete* factor- 
izability. Also, the distribution matrices follow directly 
from the discrete forms for M presented earlier. The 1 
same approach can be* applied for incompressible-flow 
problems and to staggered-grid discretizations as we'll. 

The starting point is the discretization (21). The* 
way proposed to improve* the discrete* full-potential op- 
erator is to change the discretization of Q h to Q h +A h . 
Them the discrete* full-potential operator is changed to 

Q h A h +Q h Q h -c 2 A 2, \ 
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whore A h = (g fc )‘ V V» = T" - (Q h Q h - 
c-A- h ), and T h is a desired approximation for the full- 
potential factor. In smooth regions, ,4'* is second-order 
small (proportional to /r), hence the overall second- 
order discretization accuracy is not compromised. The 

operator j is a nonlocal operator and its intro- 
duction can be effected through a new auxiliary vari- 
able v h and a new discrete equation Q h t ,h = V h p h . 
Thus, the corrected discrete approximation to (8) is 
as 

Qi> () 0 0 ydf' 0 

0 Q h 0 0 jOf 0 

0 0 Q" 0 0 

0 0 0 Q h -P' 1 0 

pr-dj h pr’df 1 pr 2 &i h 1 Q" 0 

dof ^d'i h 0 0 Q h 

(23) 

The corresponding distribution matrix. M h . for dis- 
tributed relaxation is defined as 


defined 


L h = 


M h 


• 1 0 0 0 -~ p df 0 " 

0 1 t) 0 -\0- y h 0 

0 0 1 0 -jOz h 0 

0 0 0 1 p'' o 

0 0 0 0 Q h 0 

0 0 0 0 0 1 . 


(24) 


so that the resulting matrix L h M h becomes lower tri- 
angular as 


L h M h 


Q h 0 0 0 0 0 

0 Q h 0 0 0 0 

0 0 Q h 0 0 o 

() 0 0 Q b 0 o 

pddf pddf pdd f 1 J* (» 

do?' ddf dd:" 0 Q h 


(25) 


Tlie scheme as defined above is valid for noncon- 
servative flows. A version to be used for distributed 
relaxation of conservative equations has also been de- 
signed." 

Concluding Remarks 

A general multigrid approach to attain TME solvers 
for large-scale CFD applications has been outlined. 
This approach focuses on fast convergence to the solu- 
tion of the differential equations, not necessarily pro- 
viding fast asymptotic convergence rates. The consid- 
ered measure of the method efficiency is convergence 
to the differential solution, i.e.. fast reduction of alge- 
braic errors below the discretization error level on each 
mesh. 

Because 1 the governing equations are a coupled set of 
nonlinear conservation equations with discontinuities 


and singularities, attainment of full efficiency requires 
that each of three error components be addressed: (1) 
high-frequency errors (2) uniformly smooth errors and 
(3) characteristic errors. These errors an 1 reduced 
through a combination of distributed relaxation and 
local relaxation at each grid. The relaxations are fol- 
lowed by an FAS correction from a coarser grid where 
the initial roarso-grid residuals are found from the 
fine-grid residuals of the conservative equations with 
(conservative) full-weighting restriction. 

The distributed relaxation procedure is designed to 
reduce errors in the smoothly varying regions of the 
domain and relies on the 1 factorizability property of 
the governing differential equations to isolate and treat 
optimally different factors arising in the detei ininant of 
the differential operator. Optimal relaxation of some 
particular factors may itself involve 1 a separate inner 
multigrid cycle over a limited subdomain. 

Local relaxation is a procedure designed to reduce 
large 1 residuals of the 1 conservative equations at dis- 
ce)nt,inuiti( l s/siiigulariti( l s/boundaries anel is applied in 
these regions as well as in general where the residu- 
als are largest. It employs a robust scheme (e.g., some 
block relaxation) applied before 1 and after a distributed 
relaxation sweep. 

The 1 particular factors arising in distributed relax- 
ation of the Euler and Navier-Stokes equations have 
beern discussed from the standpoint of the differential 
anel the eliscrete equations. Methods for relaxing and 
analyzing these factors within the niultigrid context. 
have 1 been presented and evaluated. Recent, progress in 
development textbook-efficient niultigrid solvers based 
ein the distributed relaxation approach has been sum- 
marized. 
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